%% Script for Generating DC Fiber Coordinates
%   Written by Bryan Howell on 1/7/2013
%   Modified and Comment on 2/25/2015
%   The trajectory of each fiber is calculated based on three reference
%   points:
%   1. vertebral level where DC fibers enter spine
%   2. vertebral level with a known topgraphic organization
%   3. vertebral level with both cuneate and gracile fasciculi
%
%   This script was used in the following PLOS ONE publication:
%   Citation: Howell B, Lad SP, Grill WM (2014) Evaluation of Intradural 
%   Stimulation Efficiency and Selectivity in a Computational Model of 
%   Spinal Cord Stimulation. PLoS ONE 9(12): e114938. 
%   doi:10.1371/journal.pone.0114938 

clear all;
clc;

%% Import Spinal Cord Geometry

patnum=5; % patient number (1-5)
cordgeom=load('geom_patientcord.txt'); % patient spinal cord geom.

% white matter boundary defined with an ellipse
ae=cordgeom(patnum,1); % transverse, left-right axis of ellipse
be=cordgeom(patnum,2); % transverse, down-up axis of ellipse
c_white=[0;cordgeom(patnum,3);0]; % center of ellipse

dr_th1=cordgeom(patnum,4);
dr_th2=cordgeom(patnum,5);
dc_th1=dr_th2;
dc_th2=cordgeom(patnum,6);

yf=cordgeom(patnum,7);

%% DC Fiber Geometry
%   assumptions:
%   - each spinal level has ~ 11 medial-lateral layers
%   - the most lateral layer contains a-beta fiber at entry + ascending
%   branch
%   - the 10 most medial layers contain DC fibers of caudal a-beta fibers

bilat=1; % 0 = unilateral, 1 = bilateral

D=9; % diameter of DC fiber (ascending branch of a-beta fiber)
nodclv=10; % number of layers w/ DC fibers
axonperlv=10; % number of DC fibers per layer
noaxons=nodclv*axonperlv; % total number of fibers/axons

%% Arc Length of DC fiber
%   NoR-MYSA-FLUT-STIN-...-STIN-FLUT-MYSA-...
%   (NoR to 2nd MYSA is 1 segment)

filename=['axonarc_',num2str(D),'um.txt'];% arc length of 1 segment
axonarc=1e-3*load(filename);% um -> mm

maxarc=280; % maximum arc length
norepeats=floor(maxarc/axonarc(end)); % number of repeats (segments)

delarc=diff(axonarc)'; % distances between successive elements in 1 segment
arcval=delarc*ones(1,norepeats); % outer product 
arcval=[0;cumsum( arcval(:) )]; % string matrix from column 1 to last column
ck=(abs(arcval)<1e-10); % values less than 1x10^-10 are considered...
arcval(ck)=0; % ...zero
noarcval=length(arcval); % total number of elements

%% 1st Reference Plane
%   - vertebral origin of 10 most medial layers
%   - e.g., see Feirabend et al. (2002)
%   - 11 layers at T11: T11 most lateral, T10-S5 split into 10 layers
%   T = thoracic
%   L = lumbar
%   S = sacral
%
%   - Assume the following:
%   1. the distance between vertebral levels is ~ 25 mm
%   2. at each vertebral level, the 11th layer is where the a-beta fibers
%   ascending branch resides

% split lateral most layer into 10 radial levels 
xyref_1=seed_dclayer(1,10,dr_th1,dr_th2,0,0.08*be,ae,be);
zref_1=-(25:25:10*25)'; % z entry point of a-beta fibers for 10 medial layers 

x1_all=xyref_1(:,1)*ones(1,10); % repeat x coord. of 10 radial lv. 10 times
y1_all=xyref_1(:,2)*ones(1,10); % repeat y coord. of 10 radial lv. 10 times
% 10 medial layers; each layer is -25 mm less than the previos one
z1_all=(zref_1*ones(1,10))'; 

% condense 3 matrices into 1 matrix
xyzref_1=cat(2,x1_all(:),y1_all(:),z1_all(:)); % (x,y,z) coord. of 1st ref.

%% 2nd Reference Plane
%   - T11

% split 10 most medial layers into 10 radial levels each
xyref_2=seed_dclayer(nodclv,axonperlv,dc_th1,dc_th2,0,yf*be,ae,be);
zref_2=zeros(noaxons,1); % xy plane at z=0

xyzref_2=cat(2,xyref_2,zref_2); % (x,y,z) coord. of 2nd ref.

%% 3rd Reference Plane
%   - C6
%   Assume:
%   1. cuneate fasciculus is present
%   => gracile fasciculus is compressed into a smaller area
%  * thus, gracile fasciculus is more medial

c6_to_t11=12; % 12 vertebral levels from C6-T11, including C6

% split CONDENSED lateral most layer into 10 radial levels 
xyref_3=seed_dclayer(nodclv,axonperlv,65,90,0.4*be,0.95*be,ae,be);
zref_3=25*c6_to_t11*ones(noaxons,1); % 25 mm/level * 12 levels

xyzref_3=cat(2,xyref_3,zref_3); % (x,y,z) coord. of 3rd ref.

%% Translation

% translate xyz coord. to center of white matter ellipse
% - this step is only crucial if coordinating with FEM model of spine
% (see PLOS ONE article)
xyzref_1=xyzref_1+( c_white*ones(1,noaxons) )';
xyzref_2=xyzref_2+( c_white*ones(1,noaxons) )';
xyzref_3=xyzref_3+( c_white*ones(1,noaxons) )';

%% Trajectory of DC Fiber
%   - based on three reference points
%   - piecewise linear: linear from ref. 1 to ref. 2 and linear from ref. 2
%   to ref. 3

dcol_coord=cell(noaxons,1); % cell array holding DC fiber population coord.

for ii=1:noaxons % for each fiber/axon

    % displacement vectors for axon at 3 ref. planes
    r1=[xyzref_1(ii,1),xyzref_1(ii,2),xyzref_1(ii,3)]; % vert. lv. of origin 
    r2=[xyzref_2(ii,1),xyzref_2(ii,2),xyzref_2(ii,3)]; % T11
    r3=[xyzref_3(ii,1),xyzref_3(ii,2),xyzref_3(ii,3)]; % C6

    r21=r1-r2; % vector that points from r2 to r1
    % update r1 so that it is half the total arc length of fiber from 2nd
    % ref. point
    r1=r2+r21/norm(r21)*(maxarc/2);
    rmag_12=norm(r2-r1); % half the distance of DC fiber

    dcol_xyz=zeros(noarcval,3); % xyz coord. for a given axon
    for jj =1:noarcval        
        if(arcval(jj)<rmag_12) % 1st half of fiber (starting from ref. 1)
            rc=r2-r1; % points from r1(new) to r2
            % travel arc length along rc
            dcol_xyz(jj,:)=r1+rc/norm(rc)*abs( arcval(jj) );
        else % 2nd half of fiber
            rc=r3-r2; % points from r2 to r3
            % travel arc length along rc
            dcol_xyz(jj,:)=r2+rc/norm(rc)*abs( arcval(jj)-rmag_12 );
        end       
    end
    dcol_coord{ii}=dcol_xyz; % place xyz coord. in cell array
    
end
 
% if bilateral, assume symmetry
if(bilat==1)
    tmp_coord1 = dcol_coord;
    tmp_coord2 = dcol_coord;
    dcol_coord = cell(2*noaxons,1);
    for k = 1:noaxons
        tmp_coord2{k}(:,1) = -tmp_coord2{k}(:,1);
        dcol_coord{k} = tmp_coord1{k};
        dcol_coord{noaxons + k} = tmp_coord2{k};
    end
    noaxons = 2*noaxons;
end

%% Storing Data in a Single Matrix

m=noaxons*noarcval; % number of elements per axon * number of axons
n=3; % three spatial dimensions
Mcoord=zeros(m,n);

% condense entire cell array into a single matrix
nonodes=(noarcval-1)/8+1;
for ii =1:noaxons
    a=(ii-1)*noarcval + 1;
    b=ii*noarcval;
    Mcoord(a:b,:)=dcol_coord{ii};    
end

Mcoord = Mcoord'; % transpose

file2 = ['p',num2str(patnum),'_dcolmxyz_',num2str(noaxons),...
         'axons_',num2str(nonodes),'nodes_',num2str(D),'um.txt'];

% % not a crucial step     
% % bounds of FEM model     
% ck1=( min( Mcoord(1,:) ) > -50 )&( max( Mcoord(1,:) ) < 50 );
% ck2=( min( Mcoord(2,:) ) > -50 )&( max( Mcoord(2,:) ) < 50 );
% ck3=( min( Mcoord(3,:) ) > -150 )&( max( Mcoord(3,:) ) < 150 );
% if(ck1&&ck2&&ck3) % if DC fiber coordinates w/in bounds of FEM model
%     dlmwrite(file2,Mcoord);
% end

%% Visualize Fibers

% unit ellipse
th_unit = 0:pi/20:2*pi;
runit   = ae*be./sqrt((be*cos(th_unit)).^2+(ae*sin(th_unit)).^2); 
xunit   = c_white(1) + runit.*cos(th_unit);
yunit   = c_white(2) + runit.*sin(th_unit);

figure;
hold on;
% reference points
plot3(xyzref_1(:,1),xyzref_1(:,2),xyzref_1(:,3),'bx');
plot3(xyzref_2(:,1),xyzref_2(:,2),xyzref_2(:,3),'bo');
plot3(xyzref_3(:,1),xyzref_3(:,2),xyzref_3(:,3),'bs');
% 2 ellipses to define 1 vertebral level
plot3(xunit,yunit,zeros(size(xunit)),'k--','LineWidth',1);
for ii = 1:length(zref_1);
    plot3(xunit,yunit,ones(size(xunit))+zref_1(ii),'k--','LineWidth',1);
end

for ii = 1:noaxons
    plot3(dcol_coord{ii}(:,1),dcol_coord{ii}(:,2),dcol_coord{ii}(:,3),'k.-');

end
plot3(dcol_coord{25}(:,1),dcol_coord{25}(:,2),dcol_coord{25}(:,3),'r.');
hold off;
title('Seeding Dorsal Column Fibers','FontSize',36,'FontWeight','b');
xlabel('x distance (mm)','FontSize',30);
ylabel('y distance (mm)','FontSize',30);
zlabel('z distance (mm)','FontSize',30);
axis([-15,15,5,35,-275,300]);
set(gca,'FontSize',30);
axis square;

% % XY Locations of DC Fibers
% xyloc = xyref_2;
% xytmp = xyloc;
% xytmp(:,1) = -xytmp(:,1);
% xyloc = cat(1,xyloc,xytmp);
% plot(xyloc(:,1),xyloc(:,2),'k.');
% dlmwrite(['p',num2str(patnum),'_dc_xyloc.txt'],xyloc);
% xyloc = xyref_2;
% xytmp = xyloc;
% xytmp(:,1) = -xytmp(:,1);
% xyloc = cat(1,xyloc,xytmp);
% [nr,nc] = size(xyloc);
% xyloc = cat(2,xyloc,zeros(nr,1));
% plot3(xyloc(:,1),xyloc(:,2),xyloc(:,3),'k.');
% dlmwrite(['p',num2str(patnum),'_dc_xyzloc.txt'],xyloc);
% return
    